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Abstract 

Photoacoustic computed tomography (PACT) is an emerging computed imaging modality that 
exploits optical contrast and ultrasonic detection principles to form images of the absorbed optical 
energy density within tissue. When the imaging system employs conventional piezoelectric ultrasonic 
transducers, the ideal photoacoustic (PA) signals are degraded by the transducers’ acousto-electric 
impulse responses (EIRs) during the measurement process. If unaccounted for, this can degrade the 
accuracy of the reconstructed image. In principle, the effect of the EIRs on the measured PA signals 
can be ameliorated via deconvolution; images can be reconstructed subsequently by application of 
a reconstruction method that assumes an idealized EIR. Alternatively, the effect of the EIR can be 
incorporated into an imaging model and implicitly compensated for during reconstruction. In either 
case, the efficacy of the correction can be limited by errors in the assumed EIRs. In this work, a joint 
optimization approach to PACT image reconstruction is proposed for mitigating errors in reconstructed 
images that are caused by use of an inaccurate EIR. The method exploits the bi-linear nature of the 
imaging model and seeks to refine the measured EIR during the process of reconstructing the sought- 
after absorbed optical energy density. Computer-simulation and experimental studies are conducted 
to investigate the numerical properties of the method and demonstrate its value for mitigating image 
distortions and enhancing the visibility of fine structures. 
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I. Introduction 

Photoacoustic computed tomography (PACT), also known as thermoacoustic or optoacoustic 
tomography, is an emerging imaging modality that holds great promise for biomedical imaging 
[?]> [?]> [?], [?], [?]• It is a hybrid modality that exploits the high optical contrast of soft tissue 
and the high spatial resolution of ultrasonic methods. In PACT, short laser pulses (typically 
nanosecond-duration) are employed to illuminate tissue. Absorption of the optical energy results 
in local heating followed by thermal expansion, which generates internal broadband photoacoustic 
(PA) wavefields via the photoacoustic effect [?], [?]. From measurements of the PA wavefields 
acquired outside the object, an image reconstruction method can be employed to estimate the 
spatially variant absorbed optical energy density within the tissue. 

When the imaging system employs piezoelectric transducers, the PA signals at the transducer 
locations are convolved with the transducers’ acousto-electric impulse responses (EIRs) during 
the measurement process. If unaccounted for, this degradation of the measurement data will result 
in a modulation of the spatial frequency components of the estimated absorbed optical energy 
density distribution [?]. In principle, the effect of the EIRs on the PA signals can be removed 
via deconvolution if the EIR is accurately known; subsequently, images can be reconstructed 
by application of a reconstruction method that neglects the EIR. Alternatively, the effect of the 
EIR can be incorporated into an imaging model and compensated for implicitly during image 
reconstruction [?], [?], [?]. 

Unlike the spatial impulse response (SIR) of a transducer, which can be described accurately 
by use of a relatively simple physics-based model [?], [?], a transducer’s EIR poses challenges 
to an analytical description [?], [?]. Various theoretical models [?], [?], [?] have been proposed 
for describing a transducer’s EIR. The parameters employed in such models, however, are either 
difficult or even impossible to measure accurately in practice. Consequently, in applications of 
PACT, it is common for the EIR to be measured experimentally. 

Although a conceptually simple task, measurement of the EIR is subject to noise and other 
errors [?], which can limit image quality in PACT. Several techniques for measuring the EIR have 
been developed [?], [?], [?], [?]. It was suggested [?] that the impulse response could be measured 
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by illuminating the transducer with an ultra-short laser pulse. However, the impulse response 
measured in this way represents the convolution of the photoacoustic pressure produced by 
parasitic sources on the surface of the transducer and the sought-after EIR [?]. Alternatively, the 
derivative of the EIR can be estimated by measuring the signal produced by optically illuminating 
an absorber that is small relative to the acoustic wavelength. In practice, signals produced by 
small absorbers can be weak [?] and errors in their low frequencies can be amplified if the 
signals are integrated to estimate the EIR. Recently, an alternate method to estimate the EIR 
was proposed to circumvent this [?], [?]. All of these methods require precise alignment of the 
acoustic source with respect to the transducer axis. When focused transducers are employed, the 
acoustic source must be aligned at the focal point. Misalignment of the acoustic source can result 
in errors in the measured EIR. In effect, the measured EIR can be contaminated by the SIR. For 
characterizing the spectral directivity of flat transducers, an optoacoustic source that produces 
quasi-plane waves was produced [?]; however, it cannot be readily utilized to characterize the 
EIR of focused transducers. Finally, when transducer arrays are purchased, although they may 
differ, the EIRs of individual elements are not typically provided, and it can be an arduous task 
to characterize each EIR. 

In this work, a joint optimization approach to PACT image reconstruction is developed for 
mitigating errors in reconstructed images that are caused by use of an inaccurate EIR. To 
accomplish this, a variable projection method [?], [?], [?] is employed to refine the measured 
EIR during the process of reconstructing the sought-after absorbed optical energy density dis¬ 
tribution. This method exploits the separable nature of the PACT imaging model. When an 
array of transducers is employed that is characterized by a collection of EIRs, the reconstruction 
method will determine a single effective EIR. Similarly, if other modeling errors are present, the 
response function produced by the method can be interpreted as an effective system response 
that minimizes the inconsistency between the measured data and the imaging model. Computer- 
simulation and experimental studies are conducted to investigate the numerical properties of the 
method and demonstrate its value for mitigating image distortions and enhancing the visibility 
of fine structures. 

The remainder of the article is organized as follows. In section II, the relevant physics and 
PACT imaging model are reviewed. The proposed image reconstruction method is described 
in Section III. The numerical studies and results are presented in Sections IV-VI. Finally, a 
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summary of the work is provided in Section VII. 

II. Background 

Below we review the basic imaging physics and discrete PACT imaging model. The reader is 
referred to [?], [?], [?], [?], [?] for comprehensive reviews of PACT. 

A. Canonical imaging model in continuous form 

In PACT, a short laser pulse is employed to irradiate an object at time t = 0 and an internal 
pressure wavefield p(r, t) is established according to the photoacoustic (PA) effect. Here, r G 
and t G [0,oo). In this work, the to-be-imaged object and surrounding medium are assumed to 
have homogeneous and lossless acoustic properties. Additionally, the width of the laser pulse 
is assumed to be negligible. Under these assumptions, the PA wavefield at a location tq G flo? 
where Qq C is the measurement aperture, satisfies 



dr. 


( 1 ) 


Here, A(r) is a compactly supported and bounded function, referred to as the object function, 
which represents the absorbed optical energy density. The quantity cq denotes the (constant) 
speed-of-sound (SOS) in the object and the background medium; /3 and Cp denote the thermal 
coefficient of volume expansion and the specific heat capacity of the medium at constant pressure, 
respectively; and V denotes the object’s support volume. 

Equation (1), which neglects the response of the imaging system as well as other physical fac¬ 
tors [?], represents an idealized imaging model for PACT in its continuous form. The associated 
image reconstruction problem is to determine an estimate of A(r) from knowledge of p{ro,t). 

B. Discrete imaging models that include transducer responses 

When piezoelectric ultrasonic transducers are employed, the photoacoustic signal p(ro,t) is 
converted to a voltage signal that is subsequently sampled. Consider the case in which the 
transducers collect data at Q locations, specified by the index g = 0,l,-- - ,Q — 1, and at each 
location S temporal samples are acquired, specified by the index s = 0,l, - - ,5' — 1. The data 
are acquired at each location with a sampling interval AT. The vector u G represents 
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a lexicographically ordered representation of the sampled voltage data, where M = QS. The 
notation [u]„ will be employed to denote the n-th element of u. 

Under the same assumptions regarding the imaging physics that are required to establish Eqn. 
(1), the measured data vector u is related to A(r) as 

[u]qs+s = Uq{t)\t=sAT = ] X / p{ro,t) dro , (2) 

jQq{rq) t=sAT 

where Uq{t) is the pre-sampled electric voltage signal corresponding to the g-th transducer whose 
active area ^q{rq) is centered at r^, and is the EIR. p(ro, t) is given by Eqn. (1). The notation 
denotes a 1-dimensional (ID) temporal convolution. Equation (2) represents a continuous-to- 
discrete (C-D) imaging model for PACT. Note that Eqn. (2) assumes that no acoustic lenses are 
attached to the piezoelectric surfaces of the transducers. 

When point-like transducers are employed, Eqn. (2) degenerates to 


[u]g5+s = h%t) *tP{rq,t) 


t=sAT^ 


( 3 ) 


where p{rq,t) is specified by Eqn. (1). 

In order to formulate the image reconstruction task as a numerical optimization problem, the 
C-D imaging model in Eqn. (2) is typically approximated in practice by a discrete-to-discrete 
(D-D) imaging model [?]. To establish a D-D imaging model, the object function A(r) can be 
approximated as 

N-l 

Aa{r) ^^[0]n(pn{r), (4) 

n=0 

where the subscript a indicates that Aa(r) is an approximation of A(r), [6]n is the n-th component 
of the coefficient vector 6, and {<pni'>^)}n=o ^^e expansion functions. In this work, interpolation- 
based expansion functions [?], [?] are employed. 

On substitution from Eqn. (4) into Eqn. (2), a D-D imaging model can be established as [?], 

[?] 

u - ne, (5) 


where the M x N matrix is commonly known as the system matrix. The system matrix H 
depends on the EIR, SIR, and the choice of expansion functions. Specifically, the elements of 
H are a function of the sampled EIR values, which will be represented by the vector h G 
Namely, [h]j = h^{iAT) for i = 0,1, • • • ,1 — 1, where / denotes the number of samples required 
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to represent the EIR. In practice, / <C S'. The explicit forms of H that is employed in this study 
are provided in Appendix A. 

To emphasize the dependence of H on h, the D-D imaging model will be expressed as 

u = H(h)6». (6) 

The accuracy of the system matrix will be degraded when the measured EIR contains errors. 
When an inaccurate system matrix is employed in an iterative image reconstruction method, 
the resulting images can contain distortions and artifacts [?]. Below, we propose a method to 
circumvent this. 

III. PACT IMAGE RECONSTRUCTION WITHOUT ACCURATE KNOWLEDGE OF TRANSDUCER 

RESPONSES 

A. Formulation of the image reconstruction problem 

We formulate image reconstruction as a numerical optimization problem 

(e,h) = argmin (/p (0, h), (7) 

6»>0,h 

where the cost function ip{6,h) is defined as 

V?(6>, h) = ||u - H(h)6l|p + \Ri{e) + «i? 2 (h). (8) 

Here, Ri{d) and i? 2 (h) represent penalty terms, whose impacts are controlled by the regulariza¬ 
tion parameters A and a, respectively. The constraint 0 > 0 in Eqn. (7) reflects that A(r) > 0 
and 0n(r) > 0 for the interpolation-based expansion functions employed in this work. If the 
expansion functions are not non-negative, this constraint should not be enforced. 

Equation (7) is fundamentally different from the conventional formulation of PACT image 
reconstruction [?], [?] in that the EIR is treated as an unknown to be estimated along with the 
approximation of A(r). This provides the opportunity for the experimentally-measured EIR to 
be refined during image reconstruction. Since Eqn. (8) is non-convex, determining the solution 
to Eqn. (7) can present challenges. As demonstrated below, the use of experimentally measured 
EIRs can provide relatively good initial estimates of h that will help the optimization algorithm 
avoid local minima. It is also important to properly design the penalties— Riid) and 7?2(h)—to 
regularize the solution. 
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B. Variable projection method 

A variable projection (VP) method [?] is employed to reformulate the minimization problem 
given in Eqn. (7). This approach is motivated by previous studies in which the VP method was 
employed successfully to estimate unknown parameters of a system matrix in separable inverse 
problems [?], [?], [?]. 

The VP method is based on the observation that 


h = argmin(/7(0,h), 


(9) 


where (0,h) is defined in Eqn. (7). Inspired by this observation, h can be parameterized as 


h*(0) = argmin(/7(0, h). 

h 

By use of this parameterization, it can be verified [?], [?] that 6 can be computed as 


( 10 ) 


0 = argmin9?(0, h*(0)), (11) 

0^0 

and, subsequently, h can be computed via Eqn. (9). In this way, the original optimization problem 
in Eqn. (7) can be solved by consideration of the two subproblems in Eqns. (9) and (11). 

It will serve useful to note that the gradient of 9?(0, h*(0)) with respect to 6 can be computed 


as 


V0p{e,h*{d)) = V0p{e,h) 


( 12 ) 


where Vg denotes the discrete gradient operator with respect to 0. The derivation of Eqn. (12) 
makes use of the optimality condition for Eqn. (10); namely, VhV^(0,h) = 0 at the point 
h = h*. Equation (10) simplifies the gradient calculation; the gradient computation prescribed 
by Eqn. (12) is identical to that employed by standard gradient descent methods for penalized 
least squares reconstruction problems. 


C. VP algorithm 

A VP algorithm for solving Eqns. (9) and (11) is provided in Algorithm 1. An experimentally- 
measured EIR, denoted by h°, is utilized to initialize h. Initialization of 9 is achieved by 
solving a constrained optimization problem in Line-2. A variety of established iterative image 
reconstruction algorithms can be employed to accomplish this [?], [?]. In Line-6, the estimate of 
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Algorithm 1 Variable projection (VP) algorithm for joint estimation of 6 and h 
1 : h° ^ experimentally-measured EIR. 

2 : 0° arg min 0 >o ^{ 0 , h°) 

3: A: 0 {A: is the number of algorithm iteration} 

4: while stopping criterion is not satisfied do 
5: k i — h 1. 

6 : arg minh h) 

7: 0 ^ 

8: end while 

9: 0 ^ 0 ^ and h ■«- h^. 


h at the k-th iteration, denoted by h^, is obtained according to Eqn. (10) by use of the previously- 
estimated 0^“^. The updating scheme for 0^ in Line-7 represents a gradient descent step for 
solving Eqn. (11), where Vc denotes the operator that projects all negative values to 0. The 
gradient is computed as specified in Eqn. (12). Note that Line-7 does not fully solve Eqn. (11). 
Instead, Line-7 moves 0^ along the negative gradient of (/? by a small step size, denoted by 7 ^. 
This distinguishes the VP algorithm from a block-coordinate descent algorithm [?], in which 
0^ argmin 0 >o (/ 2 ( 0 , h^). Note that Line-7 can be computed much more efficiently than the 
problem argmin 0 >o h^)- In addition, VP algorithms have been reported to possess faster 
convergence rates and may be less likely to be become trapped by local minima as compared 
to block-coordinate descent algorithms [?]. 

D. Implementation of the VP algorithm 

Numerical details regarding the solution of the sub-problems defined by Lines-2 and -6 in 
Algorithm 1 are provided below. 

A method for solving the constrained minimization problem in Line-2 has been described in 
[?], [?]. In this study, we assume that Ri{0) is differentiable and therefore, the constrained opti¬ 
mization problem can be solved by use of a projected gradient descent algorithm. In particular, 
we employ the updating scheme 

6 | 0 d ^ _ y V 0 <^( 6 >°’^’“\li°)}, (13) 
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where 0 °’^ denotes the estimate of 0 ° after the j-th iteration and 7 ^ denotes an updating step size. 
The step size is determined by use of a line search method [?], and the gradient is calculated as 

V0^(0,h) = {H{h)y{u-H{h)0)+VgRii0), (14) 

where (•)^ denotes the matrix transpose operator. For Ri{0) with a typical quadratic form, the 
computation of V6»i?i(0) is straightforward. Note that Eqn. (13) is of the same form as the 
updating scheme in Line-7 of Algorithm 1, suggesting the same numerical procedure can be 
employed to implement both lines. 


The second sub-problem in Line -6 of Algorithm 1 

can be efficiently implemented due to 

relatively low dimension of h—less than 100, typically. To be specific, we assume that i? 2 (h) = 

||Dh p, where D G is given by 






(l 

0 

0 ••• 

0 o\ 



-1 

1 

0 ••• 

0 0 


D = 

0 

-1 

1 • • • 

0 0 

(15) 


U 

0 

0 ••• 

-1 V 


In this case, Line -6 can be implemented 

as 




: 

= (P^P + aD^D)- 

-^P^u, 

(16) 


where the matrix P satisfies 

P( 6 »^-^)h = H(h)0'=-\ (17) 

The matrix P is described in Appendix-B. Because of its small size, P*P-|-q;D*D can be stored 
in random access memory and efficiently inverted by use of established algorithms. In the studies 
below, this was accomplished by use of the LU decomposition method [?]. 

IV. Description of computer-simulation studies 

Computer-simulation studies were conducted to investigate the numerical properties of the VP 
algorithm. 
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A. Simulation of noise-free data 

The numerical phantom shown in Figure la was employed. The phantom had a support area of 
22.0 X 22.0 mm^ and contained six uniform disks that were assigned different values of absorbed 
optical energy density. 

A 2D circular measurement geometry was employed. Q = 128 transducers were evenly 
distributed on a ring of radius 25 mm that enclosed the phantom. The SOS was assumed to 
be constant and set at cq = 1.5 mm//is. Since the simulated data were formed by use of the 
C-D imaging model in Eqn. (2), no inverse crime was committed. The components of this 
vector corresponded to T = 600 equally spaced temporal samples over the interval [10,25) /us. 
Subsequently, the noiseless voltage vector Ug was obtained by convolving the pressure data with 
EIR-1 in Figure lb. 

B. Simulation of noisy data 

Noisy measurement data were computed as 

u Ug + ri, (18) 

where ri is a random vector whose components were independent and identically distributed 
Gaussian random variables. The standard deviation of each element of n was 3% of [ugjmax! 
where [ugjmax denotes the maximum value contained in Ug. 

C. Implementation of image reconstruction algorithms 

From the simulated noiseless and noisy data, images were reconstructed by solving the 
minimization problem in Eqn. (7) by use of Algorithm 1. Conventional quadratic smoothness 
penalties were employed: 

N-l 

n=0 keAfn 
1-2 

i?2(h) = 5]([h],+i - [h]if, (20) 

i=0 

where Mn is the index set of four neighboring pixels of the n-th pixel. 

The reconstruction region (22.0 x 22.0 mm^) was represented by 440 x 440 pixels with pixel 
size 0.05 mm in each dimension. The initial guess of the EIR employed in the VP algorithm 
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was different than the EIR that was assumed when generating the simulated data. This served 
to simulate a situation in which an experimentally measured EIR contained errors. 

Each element in a real-world transducer array possesses its own EIR. In practice, the dif¬ 
ferences between the EIRs are sometimes neglected and an EIR corresponding to a single 
element may be used to represent all elements in the array. In some of the studies below, 
the EIR employed to initialize the VP algorithm (EIR-2 in Figure lb) and the EIR employed 
to produce the simulated measurements (EIR-1 in Figure lb) were experimentally measured 
from two different transducer elements in a circular transducer array (see Sec. VI-B). EIR-1 was 
measured by temporally integrating the PA signal produced by a point source positioned at the 
focus of the transducer. EIR-2 was measured by use of the method reported in [?]. In order 
to investigate the sensitivity of the VP algorithm to the initialization of the EIR, we employed 
different EIRs obtained by degrading EIR-1 as described later. When solving the sub-problem 
in Line-2 of Algorithm 1, 6 was initialized as the zero vector. Algorithm 1 was terminated after 
500 iterations, since it was observed that the changes in the reconstructed images with more 
iterations were negligible. When implemented by use of a single core of an Intel Xeon E5-2640 
CPU, each iteration required approximately 7s to complete. 

For comparison, we also reconstructed images by use of a conventional gradient-based iterative 
image reconstruction algorithm that considered the EIR to be fixed. This algorithm was the same 
as the one employed to compute the initial guess of 6 in Line-2 of Algorithm 1, which was 
described in Section III-D. As with Algorithm 1, each iteration required approximately 7s to 
complete. The reconstruction algorithm was run for 150 iterations, since the changes in the 
reconstructed images with more iterations were negligible. Note that the computational cost of 
Line-6 in Algorithm 1 was about 5% of the Line-7 in Algorithm 1, which is why each iteration 
took almost the same time in both the conventional iterative method and the VP algorithm. 

D. Image accuracy assessment 

The accuracy of the reconstructed images was assessed in terms of the root-mean-squared-error 
(RMSE) between the reconstructed image and the true phantom as 



( 21 ) 
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where N is the number of pixels, and [9]n and [0o]n the n-th pixel values of the reconstructed 
image and phantom, respectively. 

V. Computer-simulation results 
A. Images reconstructed from noise-free data 

1) Mitigation of artifacts and distortions caused by errors in the assumed EIR: Figure 2a 
shows the image reconstructed by use of the conventional iterative method that utilized a system 
matrix based on EIR-2. Different values of the regularization parameter A from the interval 
[0,10“^] were considered. The reconstructed image with the value of A that minimized the RMSE 
was chosen to represent the best performance of the conventional iterative method. Eigure 2a 
and the profile in Eigure 2c demonstrate that the use of an inaccurate EIR can result in strong 
artifacts and distortions in images reconstructed by use of the conventional methods. 

When the VP algorithm was applied, different values of the regularization parameter A from 
the interval [10“®,10“^] and a from the interval [200,20000] were considered. The image that 
minimized the RMSE was chosen and displayed in Eigure 2b. As revealed by this image and 
the profiles in 2c, the VP algorithm yielded an image with fewer artifacts and distortions, and 
image fidelity was improved as reflected by the reduced RMSE. 

2 ) Effect of frequency contents of the objects and EIR: Since the voltage signal is generated 
through convolution of the pressure data and EIR, the EIR serves as a bandpass filter. Thus, the 
information contained in the high frequency components is lost in the resulting voltage signal. 
We conducted a series of computer-simulations to show that the accuracy of the reconstructed 
9 and h will be affected by this loss of information. 

The original sharp phantom shown in Eigure la was convolved with a Gaussian blurring 
kernel to generate a smoothed phantom that possessed smaller relative spatial bandwidths. We 
employed the sharp and the smoothed phantoms to generate pressure data; the pressure data 
generated by the sharp phantom had a larger bandwidth than that generated by the smoothed 
one, as shown in Figure 3. 

The results shown in Figure 4 suggest that the reconstructed estimates of the EIR become 
more accurate when the bandwidth of the A(r) is increased (Figure 4g and 4h). On the other 
hand, the reconstructed estimates of A(r) become more accurate when the bandwidth of the 
EIR is increased (Eigure 4a and 4b). Eor a given EIR, the reconstructed estimates of A(r) that 


May 28, 2015 


DRAFT 


13 


contain sharp features contain more oscillations than the estimates corresponding to the smoothed 
versions of ^(r). This is because more high frequency information is lost during the convolution. 

3) Effect of data incompleteness: Incomplete, or sparsely sampled, data sets are sometimes 
acquired in practice. To study the effect of data incompleteness on the VP algorithm, we 
reconstructed images from data corresponding to half of the equally spaced transducers (Q = 64). 
Because the data were noiseless, no explicit regularization was employed (A = 0) in the 
conventional reconstruction algorithm. However, the explicit regularization was still employed 
in the VP algorithm because of the ill-posed nature of the joint reconstruction problem. The 
results are shown in Figure 5. As expected, use of the incomplete data set resulted in less 
accurate reconstructed images for both the conventional iterative reconstruction method and the 
VP algorithm. However, this effect was more pronounced for the VP algorithm. Note that for the 
VP algorithm, larger values of the regularization parameters were applied when the incomplete 
data set was employed than when the complete data set was employed (Figure 5h and 5g). 

4) Effect of initial estimate of EIR: The robustness of the VP algorithm with respect to 
perturbations in the EIR was investigated. Perturbed EIRs were generated by adding different 
levels of random noise to the low frequency components (first 10% of the total bandwidth, except 
for the DC component) of the true EIR. The similarity of a perturbed EIR to the true EIR was 
quantified by the correlation coefficient, which is defined by 

(hi - /ihi)^(h2 - yUhJ 


P = 




( 22 ) 


where cTh^, is the standard deviation of h^, /ih^ is the mean of h^, k = 1,2, and / is the length 
of hi and h 2 . The value of p ranges from —1 to 1. The maximum value of p is achieved when 
one EIR is linear with respect to the other EIR with a positive slope (i.e., hi = ah 2 + b, for 
some constant a > 0 and 6), which indicates that the two EIRs are ‘identical’ to each other in 
terms of similarity. On the other hand, p equals —1 when hi = ah 2 + 6, for some constant a < 0 
and b. 

As shown in Figure 6, when the error in the EIR was small (e.g., as with the EIR in Figure 6a), 
images were reconstructed with high accuracy using the VP algorithm. When the perturbations 
in the EIR were stronger (e.g, as in Eigure 6c), artifacts and distortions in the reconstructed 
images were still significantly reduced by use of the VP algorithm; however, larger values of 
the regularization parameters had to be applied. When p < 0 as in the initial EIR in Figure 6e, 
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no improvement was observed in the image reconstructed by use of the VP algorithm. 

B. Images reconstructed from noisy data 

1) Mitigation of artifacts and distortions caused by an inaccurate EIR: Figure 7a reveals 
that use of the inaccurate EIR in the conventional iterative method created strong artifacts and 
distortions. Figures 7b confirms that the artifacts and distortions were significantly mitigated 
when the VP method was employed. Image profiles for both cases are shown in Figures 7c. 
The overall accuracy of the recovered EIR, shown in Eigure 7d and It, was improved, but it 
contained spurious oscillations. 

2 } Continuous dependency on regularization parameters: Images reconstructed by use of the 
VP algorithm with different values of the regularization parameter values are shown in Figure 8. 
The recovered EIRs and their corresponding Fourier spectra are shown in Figures 9 and 10, 
respectively. The RMSE values are computed and displayed together with the corresponding 
images. As expected, the images reconstructed with smaller values of A contain higher noise 
levels, while images using larger A possess a reduced noise level. However, larger values of 
A also caused artifacts in the reconstructed images. The same observation can be made for 
the effect of the regularization parameter a on the recovered EIR. One also observes that the 
reconstructed images and EIRs depend continuously on the regularization parameters A and a, 
i.e. small changes in the regularization parameters cause minor changes in the reconstructed 
images and EIRs. 


VI. Experimental validation 

The proposed algorithm was further investigated by use of experimental data acquired from 
a 512-element full-ring-array photoacoustic computed tomography system [?]. 

A. Phantom objects 

The first phantom was comprised of a single black needle of diameter 0.25 mm and length 20 
mm embedded in an agar gel. The second phantom was comprised of a mouse kidney embedded 
in an agar gel. In both experiments, the phantom and the transducer array were aligned so that 
the object of interest laid in the focal plane of the transducer array. 
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B. Data acquisition 

The illumination light source was a tunable optical parametric oscillator (OPO) laser (basiScan, 
Spectral Physics) pumped by a Nd:YAG laser (Brilliant B, Quantel) with 12 ns pulse duration and 
10 Hz pulse repetition rate. Before reaching the sample surface, the laser beam was homogenized 
using an optical diffuser (EDC-5 Photonics) to form a 25-mm-diameter circular light beam. For 
both experiments, we tuned the OPO laser output to a wavelength of 610 nm . The maximum 
light intensity at the surface of the sample was approximately 5 mJ/cm^, which is well below 
the American National Standard Institute (ANSI) limit at 610 nm (20 mJ/cm^) [?]. 

The PA signals were detected by a 512-element full-ring transducer array with 5 MHz central 
frequency (80% bandwidth) and 50 mm ring diameter. Each element in the array was mechani¬ 
cally shaped into an arc to produce an axial focal length of 19 mm. At each transducer location, 
1300 temporal samples were acquired at a sampling rate of 40 MHz. Accordingly, the dimension 
of the measured data set was 1300 x 512. Additional details of the imaging system can be found 
in [?], [?]. 

C. Image reconstruction 

Two numerical imaging models were employed in the studies involving experimental data. 
Both models are described in Appendix A. In the first, the SIR effect was not considered and a 
2D interpolation-based D-D imaging model was employed. Most of the presented results were 
reconstructed by use of this imaging model. For the needle phantom, the size of the reconstructed 
region was 22.0 x 22.0 mm^, which was represented by 440 x 440 pixels. For the mouse kidney, 
the size of the reconstructed region was 14.0 x 14.0 mm^, which was represented by 280 x 280 
pixels. For both objects, the pixel size was 0.05 x 0.05 mm^. The needle and kidney data were 
processed using a single core of an Intel Core 17-3770 CPU (4 cores, 3.4 GHz). It required 
approximately 10.5 s and 5 s to complete one iteration for the needle and kidney data sets, 
respectively. The second imaging model included SIR effects and was based on a 3D spherical 
voxel imaging model [?], [?]. This model was only applied to the kidney data. In this case, only 
a single slice through the object at the focal plane was reconstructed. The object was represented 
by 280 X 280 x 1 voxels of dimension 0.05 x 0.05 x 0.05 mm^. In order to model SIR effects, 
each transducer surface was divided into four equal parts in the elevation direction, and a far- 
field approximation was employed to calculated the measured pressure for each sub-element. 
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also referred to as a ’patch’ [?]. We selected 4 patches because the maximum phase error in 
this case was less than one-eighth of the wavelength corresponding to the central frequency 
of the transducer. The algorithm was implemented by use of the CUDA parallel programming 
framework and executed on a single GPU (Tesla K20c). Processing the kidney data required 15 
minutes per iteration. Details on how the 3D model was constructed can be found in Appendix 
A. 

In both cases, images were reconstructed by use of both the VP algorithm and the conventional 
algorithm described previously. The VP algorithm was terminated after 120 iterations, while the 
conventional method was terminated after 50 iterations. The initial guess for the EIR was an 
experimentally-measured EIR from an element in the PACT system, and the initial guess for 
6 was all zeros. The regularization functions employed corresponded to those in Eqns. (19) 
and (20). The values of the regularization parameters A and a were determined empirically. 
We swept the values of these parameters over wide ranges with a small step size. Instead of 
attempting to identify optimal regularization parameter values, which are application dependent, 
we investigated how the regularization parameter values affect the reconstructed images. 

D. Results: Needle phantom 

Figure 11 displays images of the needle phantom reconstructed by use of the simple backpro- 
jection method [?]. Figures 12 and 13 display the images reconstructed by use of the conventional 
iterative method and VP algorithm, respectively. 

Figures 12 and 13 show that the width of the needle in the reconstructed image increases as 
the regularization parameter A increases for both the conventional iterative method and the VP 
algorithm. The images reconstructed by use of the VP algorithm appear to have a reduced noise 
level compared to the images reconstructed by the backprojection and conventional iterative 
methods, regardless of the choice of the regularization parameter values. The profile plots 
corresponding to these three methods are shown in Figure 14. Since the image of the coefficient 
vector 6 and the EIR h are recoverable only up to a multiplicative constant, every profile was 
normalized for comparison. These plots demonstrate that the image reconstructed by use of the 
VP algorithm possessed a more uniform background than those obtained by the backprojection 
and the conventional iterative methods. 
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E. Results: kidney phantom 

The images and EIRs reconstnicted by use of the VP algorithm that was based on the 2D 
imaging model that neglected the SIR are shown in Figures 15 and 16. The latter figure contains 
results corresponding to different values for the regularization parameter A. From Fig. 15, it can 
be observed that use of the conventional iterative method that utilized the measured EIR resulted 
in distortions and loss of details in the reconstructed images. Use of the VP algorithm improved 
the contrast and the details in the reconstructed images (Fig. 15c and 16a). Furthermore, the 
images reconstructed by use of the VP algorithm had a more uniform background. 

In Figure 17, the results corresponding to use of the 3D imaging model that incorporated SIR 
effects are shown. The EIR estimated by the VP algorithm is also shown. In Figure 18, images 
and EIRs reconstructed by use of the VP algorithm with different regularization parameters 
values are shown. Similar to the case described above where the transducer SIR was neglected, 
these results reveal that use of the VP algorithm can produce images with a cleaner background 
and enhanced spatial resolution than yielded by use of a conventional iterative algorithm that 
employed the measured EIR. For example, detailed information regarding the vessels near the 
organ’s periphery was better preserved by the VP algorithm than by the conventional iterative 
algorithm. These images corroborate our assertion that the VP algorithm can significantly reduce 
the artifacts and distortions in the reconstructed image. It is also worth pointing out that, unlike 
the numerical phantom studies, the artifacts and distortions in the images may be caused not 
only by the inaccurate EIR but also by other factors, such as neglecting acoustic heterogeneities 
and the variation of the EIRs among the elements of the transducer array. In such cases, the EIR 
estimated by the VP algorithm represents an effective system impulse response that minimizes 
the inconsistency between the measured data and the imaging model. 

F. Auto-focus capabilities 

Conventional PACT reconstruction algorithms assume that the medium is described by a 
constant speed-of-sound (SOS) value. In practice, this value may not be known precisely and 
can be tuned [?] to maximize the spatial resolution of the reconstructed images. The effect 
of an incorrect SOS value can sometimes be compensated for by use of the VP algorithm 
due to modification of the EIR during the joint estimation. Figures 19a and 19b show images 
reconstructed by use of the conventional iterative method and the VP algorithm, respectively. 
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when different constant SOS values are assumed. The 2D imaging model that ignored the SIR was 
employed. Nearly identical images were reconstructed by use of the VP algorithm, even though 
the assumed SOS values were different in each case. The images contained reduced artifact 
levels as compared to those reconstructed by use of the conventional method. The recovered 
EIRs differed by a time shift (as displayed in Fig. 19c). Since the object was located near the 
center of the transducer array and was small compared to the radius of the array, the scaling 
effect due to the inaccurate SOS can be approximated by the shift of the EIR, which explains 
how the recovered EIR compensates for the error in SOS value. 

VII. Conclusions and discussion 

In this study, we proposed a joint reconstruction approach for PACT that mitigates artifacts 
in the reconstructed images caused by use of an inaccurate EIR. A nonlinear least squares 
minimization problem was formulated, which exploited the bi-linear structure of the imaging 
model, and a VP algorithm was employed to solve the minimization problem. The numerical 
properties of the VP algorithm were also investigated. The results demonstrate that the joint 
reconstruction approach for estimating both the system response and the absorbed optical energy 
density can increase the fidelity of the reconstructed image. Although not presented, we also 
conducted computer-simulation studies based upon an existing three-dimensional small animal 
imaging system [?], [?], and the results were consistent with those presented. 

It should be emphasized that the recovered EIR, in general, is not equivalent to the actual 
EIR of a system. Instead, the VP algorithm finds the linear temporal filter that best matches 
the measured pressure to the modeled pressure in a penalized least squares sense. If the EIR 
is the only source of model error, the filter will correspond to the EIR. However, if other 
system inconsistencies, such as sound speed variations, acoustic absorption, or the spatial impulse 
responses of the transducers, are present, the VP algorithm will produce an estimated filter that 
attempts to mitigate these sources of model error. In practice, it can be difficult or overly time- 
consuming to explicitly account for all these potential sources of inconsistency in a PACT 
reconstruction algorithm. Further, including them can result in a tremendous increase in the 
computational cost of the algorithm. Since the VP algorithm can provide a rough correction for 
these effects, it can serve as a cheap and effective way to compensate for model mismatch. 
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The minimization problem defined in Eqn. (7) is non-convex. Hence, the optimization al¬ 
gorithm may converge to a local minimum. However, the literature [?] suggests that the VP 
algorithm is more likely to converge to the global minimum than other algorithms such as block- 
coordinate descent algorithms. Our computer-simulation studies revealed that the VP algorithm 
consistently converged to accurate solutions, suggesting that utilizing proper regularization meth¬ 
ods and good initial guesses will improve the ability of the algorithm to avoid local minima. 
The experimental results confirmed the effectiveness of the proposed algorithm for mitigating 
image artifacts and distortions. 

There remain several topics for future investigation. Our current implementation involves two 
regularization parameters. Although numerical methods—such as the L-curve method [?] — 
have been proposed for determining reasonable values for these parameters, these methods do 
not work perfectly in all applications. In this study, to reveal the impact of parameter settings on 
reconstruction algorithm performance, we reconstructed a collection of images using different 
regularization parameter values. The optimal regularization parameter values should depend on 
a specified diagnostic task and observer [?] and their determination represents a topic for future 
investigation. 

There also remains a need to investigate methods for incorporating additional a priori infor¬ 
mation regarding the EIR into the reconstruction problem. If we assume the EIR is sufficiently 
smooth, spline functions are a natural choice to parameterize the EIR and reduce the number 
of unknowns in the minimization problem. We conducted numerical studies to evaluate this. 
Although not shown here, the results suggest that, depending on the interpolation points, the 
number of unknowns employed to represent the EIR can be reduced (from 64 to 32 in this 
study). However, the reconstructed images and EIRs were similar to the results obtained without 
using the spline functions. Besides, no computational advantages (such as time and memory 
usage) were observed. It is also possible to employ an analytic parameter-based EIR model [?]. 
To accurately model a realistic transducer, tens of parameters are needed. How to effectively 
solve the associated minimization problem remain s a topic for future work. Non-smooth sparsity- 
promoting penalties, such as TV penalties [?], can be applied to the absorbed energy density 
[?]. In the VP algorithm, the updating scheme for 6 is based on a gradient-descent method 
that exploits the differentiability of the smoothness penalty (i.e. Eqn. (19)). When non-smooth 
penalties are adopted, this gradient-descent method can potentially be replaced by a proximal 
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gradient algorithm [?]. 
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Appendix A 

Explicit forms of system matrices 
A. System matrix based on interpolation expansion functions 

In the 2D computer-simulation studies, an interpolation-based image model was employed. In 
interpolation-based D-D imaging model the coefficient vector is defined as samples of the object 
function on the nodes of a uniform Cartesian grid: 


[dr 5{r - rn)A{r) 


n = 0,l,--- ,N -1 


(23) 


Jv 


where r„ = (a:„, UnY specifies the location of the n-th node of the uniform Cartesian grid. The 
definition of the expansion function depends on the choice of interpolation method. If a trilinear 
interpolation method is employed, the expansion function can be expressed as 



(1 - ^A7^)(1 - \f\x - Xnl\y - yn\ < 


(24) 


otherwise 


0 , 


where is the distance between two neighboring grid points. 

In principle, the interpolation-based D-D imaging model can be constructed by substituting 
Eqns. (23) and (24) into Eqn. (2). In practice, however, implementation of the surface integral 
over Sq is difficult for the choice of expansion functions in Eqn. (24). Therefore, utilization of 
the interpolation-based D-D model commonly assumes the transducers to be point-like. 

Since Eqn. (1) can be reformulated as the well-known spherical Radon transform (SRT) 

5 f(ro,t)= [drA{r)S{cot-\ro-r\), (25) 

Jv 

where the function g{rQ,t) is related to p(ro,t) as 



(26) 


the implementation of H is decomposed as a three-step operation: 


u = H6> = H'^DG6> 


(27) 


May 28, 2015 


DRAFT 






21 


where G, D, and H® are discrete approximations of the SRT (Eqn. (25)), the differential operator 
(Eqn.(26)), and the operator that implements a temporal convolution with EIR, respectively. G 
was implemented in a way that is similar to the ‘ray-driven’ implementation of Radon transform 
in X-ray CT, i.e., for each data sample, we accumulated the contributions from the voxels that 
resided on the spherical shell specified by the data sample. By use of Eqns. (4), (23), (25), and 
(24), one obtains 


[GS] 


qS+s 




w 


qS+s^ 


(28) 


N-l Ni-lNj-l 

n=0 i=0 j=0 

where [gj^s+s ~ fi'(rg,f)|t=fcAt with specifying the location of the g-th point-like transducer, 
and Ni and Nj denote the numbers of divisions over the two angular coordinates of a local 
spherical coordinate system.The differential operator in Eqn. (26) is approximated as 


[Dg], 


f[s]gK+k+l [g]giV+fc-l 


-) - [p], 


(29) 


\qK+k ~ gTrCpAf V k + 1 k-1 J ~ ^^iqK+k' 

where [pj^s+s ~ p(rg,f)|t=sAf Finally, the continuous temporal convolution is approximated by 
a discrete linear convolution as 


[H‘P],S+, 

where [h], = 


S-1 

^^[h]s-l-ft[p]q5-|-«; = [u]q5-|_s, 
K=0 


(30) 


B. 3D spherical voxel-based imaging model including SIR 


Eor the 3D spherical voxel-based model, the expansion functions were defined as 


{ 1, if|r-r„|<e, 

0, otherwise. 


(31) 


where r = {xn,yn,Zn)^ specifies the coordinate of the n-th grid point of a u nif orm Cartesian 
lattice, (•)^ denotes the transpose of a vector, and e is the half spacing between lattice points. 
The coefficient vector 6 es defined as 


[0]n = -^^ [ (pnir)A{r)dr, (32) 

^voxel Jv 

where Vcube and V^oxei are the volumes of a cubic voxel of dimension 2e and 0n(i'), respectively. 
Let M“(f) denote the pre-sampled voltage signal that would be produced by Aa{r), where Aa{v) 
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is the approximation of A(r) established by use of the chosen expansion functions. By use of 
Eqns. (2), (4), and (31), it can be verified that 

1 ^-1 

Uq{t) = *tPo{t) n '^[9]nhl{rn,t). 


(33) 


1 n=0 


Here, po{t) is the ‘N’-shaped profile produced by a uniform sphere of radius e: 


Po{t) = 




P L 


+ -Hlt-- 
Coj V Co 


(34) 


where H{t) is the Heaviside step function and 


d{t — 


|r'-rnh 


CO 


27r|r' 


dr' 


(35) 


is the SIR of the g-th transducer. By temporally sampling (33) and employing the approximation 


u 


qS+s 




=sAT' 


the D-D imaging model in Eqn. (5) is established [?], where 


[H]g5+3,n = nPo{t) ^h^(rn,f) 


(36) 


t=sAT 


Here, [■]m,n denote the entry in the m-th row and n-th column of the matrix. When the transducer 
has a flat and rectangular detecting surface of area a x b, under the far-fleld assumption, the 
temporal Fourier transform of the SIR is given by [?] 


ki^nP) 


o6exp(-j27r/^^^2^ 

27r|r' - r„| 


Sine 


ax: 


tr 

nq 


co|r^ - 


Sine 


u^nq 


co|r^ - 


(37) 


where and specify the transverse coordinates in a local coordinated system that is centered 


nq 


about the gth transducer. 

Since the surfaces of the focused transducers employed in the reported experimental studies 
are curved, direct use of the far-field approximation assuming a flat transducer can result in 
patterned image artifacts. To alleviate this limitation, we adopt a simple divide-and-integrate 
algorithm [?], where each transducer element face is divided into m x 1 identical patches. Each 
patch is considered to be flat and described by the far-field approximation. Let ^^ j(r„,f) be 
the resulting SIRs that are specified by the patch index z = 1, • • • , m. The SIR for the original 
transducer face hg{rn,t) is then approximated by averaging the patch SIRs over all patches: 


^ m 

i=l 


(38) 
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Appendix B 

An equivalent reformulation of the imaging model 


First observe that a D-D model without considering EIR can be derived as 


Po 

Pi 


= Hp0. 


(39) 


L pq-1 J 

Here, the vector p G represents a lexicographically ordered representation of the sampled 
pressure data, the dimension P is defined by the product of the number of pressure temporal 
samples (T) acquired at each transducer location and the number of transducer locations (Q), 
and Pg G E^ (g = 0, 1 , • • • , Q — 1 ) is the sampled pressure data corresponding to location index 
q. The system matrix Hp, without considering EIR, is of dimension P x N, whose elements are 
defined by Eqn. (28) and (29). To update h using (10), another equivalent formulation of the 
D-D image model ( 6 ) can be established as 



Po 


Po n h 


Poh 

u = Hp 0 *th = 

Pi 

*t h := 

Pi n h 


Pih 


. PQ -1 . 


Pq_i h 


1 

1 


(40) 


Po 

Pi 


h = Ph. 


[Pq-iJ 


(41) 


Here, *t denotes the discrete temporal convolution and Pi is the convolution matrix corresponding 
to Pi. Matrix P is defined by Eqn. (41). The number of temporal samples S (of each transducer), 
/ (of EIR), and T (of the pressure) satisfy the relation S = I+ T — 1. With this reformulation, 
one has f{d, h) = ||u — Ph|p. 
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(a) Phantom 


—Measured EIR 1 
---Measured EIR 2 



■0.5' ^ 

0 0.4 0.8 1.2 1.6 

Time (|li s) 


(b) EIRs 


Fig. 1: (a) The absorbed energy density map employed in the computer-simulation studies, (b) 
EIRs employed for 2D experiments. 
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Fig. 2: (a) and (b) Images reconstructed from noiseless data using the conventional iterative 
method with A = 1.0 x 10“^ and the VP algorithm with A = 1.0 x 10“^ and a = 2000, 
respectively, (c) Image profiles through noiseless images seen in Fig. 2 corresponding to the 
reconstructions with the VP algorithm (solid line), with the conventional iterative method (dashed 
line), and phantom (dash-dot line). The locations of the profiles are indicated by the “Y” arrows 
in the Fig. 2a, and 2b, respectively. 
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Fig. 3: Spectrums of the pressure data generated by the sharp and smoothed objects. 
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Fig. 4: The first row shows the reconstructed image and the EIR from the smoothed object, 
where the spectrums of the generated pressure data are smaller than that of the EIR. The second 
row shows the ones from the sharp object, where the spectrums of the generated pressure data 
are larger than that of the EIR. Fig. (a, e) are the reconstructed images. Fig. (b, f) are the profile 
plots. Fig. (c, g) are the recovered EIRs, and Fig. (d, h) are the recovered EIRs in the frequency 
domain. 


May 28, 2015 


DRAFT 












FIGURES 


28 




(b) A = 0 



— Half data 


---Full data 


[A 

- - Phantom 

-i 

.n 


-11 0 11 
Y(mm) 

(c) 






Y(mm) 

(i) 


Fig. 5: Images reconstructed from incomplete data sets. The first row shows the images 
reconstructed using the true EIR and the conventional iterative method. The second row shows 
the images reconstructed using the wrong EIR and the conventional iterative method. The third 
row shows the images reconstructed using the VP algorithm. Images from left to right in each 
row are reconstructed using incomplete (half) data, full data, and the profile plots, respectively. 
The locations of the profiles are indicated by the “Y” arrows in the images. 


May 28, 2015 


DRAFT 










FIGURES 


29 




(c) p = 0.7271 


(d) A = 2.0 X 10-^Q^ = 2000 



(e) p = -0.3177 


(f) A = 10-^, a = 2000 


Fig. 6: Images reconstructed by use of different initial guesses of EIR. Each row shows the 
results using the same initial EIR. The first column shows the true , initial , and recovered EIR. 
The second column shows the images reconstructed using the VP algorithm. 
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Time (ji s) Frequency (MHz) 

(d) (e) 

Fig. 7: (a) and (b) Images reconstructed from noisy data using the conventional iterative method 
and using the VP algorithm, respectively, (c) Image profiles through the noisy images in Fig. 7a 
and 7b, corresponding to the reconstructions using the VP algorithm (solid line), the conventional 
iterative method (dash-dot line), and the phantom (dashed line). The locations of the profiles 
are indicated by the “Y” arrows in the Fig. 7a, and 7b, respectively, (d) The recovered EIR 
corresponding to the reconstruction using noisy data with the VP algorithm (solid line), the true 
EIR (dashed line), and the initial EIR (dash-dot line), (e) The spectrums of the corresponding 
EIRs. 
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(g) A = 10-^ a = 2.0 X IQS, (h) A = a = 2.0 x lO^, (i) A = IQ-^, a = 2.0 x 10^ 

RMSE = 0.0807 RMSE = 0.0353 RAISE = 0.0270 


Fig. 8: Images reconstructed from noisy data by use of the VP algorithm corresponding to 
regularization parameters A = {10“^, 10“^, 10“^} and a = {2.0 x 10^, 2.0 x 10^, 2.0 x 10^}. 
RMSE values are also shown. Images are displayed in their full dynamic ranges respectively. 
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(a) A = 10“®, a = 2000 


(b) A = 10-3, a = 2000 


(c) A = 10-\ a = 2000 




(d) A = 10-3, a = 20000 


(e) A = 10-3, a = 20000 


(f) A = 10-1, a = 20000 




0.5 


-0.5 



—Recovered 
True EIR 
Initial EIR 





Time (in s) 


(g) A = 10-3, a = 2.0 X 103 


(h) A = 10-3, a = 2.0 X 103 


(i) A = 10-1, a = 2.0 X 103 


Fig. 9: The recovered EIRs corresponding to the reconstructed images in Fig. 8. A = {10 
10-^ 10-^} and a = (2.0 x 103,2.0 x 10^2.0 x 10^}. 
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(a) A = 10-^ a = 2000 




Frequency (MHz) 


Frequency (MHz) 


(d) A = 10“®, a = 20000 



(e) A = 10-3, a = 20000 




Frequency (MHz) 


(g) A = 10-®, a = 2.0 X 10® 


(h) A = 10-3, a = 2.0 X 10® 


(i) A = 10-1, a = 2.0 X 10® 


Fig. 10: The recovered EIRs in the frequency domain corresponding to the reconstructed images 
in Fig. 8. A = {10-^ 10-^ IQ-^} and a = (2.0 x 10^, 2.0 x 10^ 2.0 x 10^}. 
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Fig. 11: Image of a needle phantom reconstructed using the backprojection. The zoomed-in 
image corresponds to the ROI of the dashed rectangle. 
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Fig. 12: Reconstructed needle phantom image using the conventional iterative method with the 
non-negativity constraint, (a) A = 1.0 x 10“^ and (b) A = 1.0 x 10“^. The zoomed-in image 
corresponds to the ROI of the dashed rectangle. 
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(b) A = 1.0 X 10"^, a = 5000 


Fig. 13: (a) and (b) Images reconstructed by use of the VP algorithm corresponding to the initial 
EIR guesses shown in the plots of Fig. 13c. (a) A = 1.0 x 10“® and (b) A = 1.0 x 10“^. The 
zoomed-in image corresponds to the ROI of the dashed rectangle. Images are displayed in their 
full dynamic ranges respectively. Fig. 13c shows the recovered and initial EIRs. 
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- 1.5 - 0.75 0 

X(mm) 


Fig. 14: Image profiles at y = 2.4 mm along the x-axis from —1.5 mm to 0 mm extracted 
from images (22.0 x 22.0 mm^ and centered at (0.0, 0.0)) reconstructed by the backprojection 
algorithm (dashed line), the conventional iterative method (CIM) with A = 1.0 x 10“^, and the 
VP algorithm with A = 1.0 x 10“^ and a = 5000. 
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(a) A = 0.0 


(b) A = 0.0 


0.3 


0.15 


0 



0.16 



(c) A = 1.0 X 10“®, a = 5000 


Fig. 15: (a) Image reconstructed by use of the conventional iterative method without the non¬ 
negativity constraint, (b) Image reconstructed by use of the conventional iterative method with 
the non-negativity constraint, (c) Image reconstructed by use of the VP algorithm, with the initial 
guess of the EIR shown in the right plots of subfigure (c). The zoomed-in image corresponds to 


the ROI of the dashed rectangle. Images are displayed in their full dynamic ranges respectively. 


The right plot shows the recovered and initial EIR. The SIR was ignored in these studies. 
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(a) A = 1.0 X 10-5, a = 5000 



— Recovered 
- Initial EIR 



Time (|j, s) 


(b) A = 1.0 X 10-4, a = 5000 



0.16 



— Recovered 
- Initial EIR 



Time (|j, s) 


(c) A = 1.0 X 10-3, a = 5000 


Fig. 16: The left panels of subfigures (a), (b), and (c) display images reconstructed by use of 
the VP algorithm and different regularization parameters with the initial guess of the EIR shown 
in the corresponding right panels. The SIR was ignored in all cases. The zoomed-in image 
corresponds to the ROI of the dashed rectangle. The grayscale windows were [0,0.17]. 
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(a) A = 1.0 X 10-5 


0.06 


0.03 


0 



(b) A = 1.0 X 10-5,0; = 5000 


Fig. 17: (a) Image reconstructed by use of the conventional iterative method with a non-negativity 
constraint. The grayscale window was [0,0.06]. (b) Image reconstructed by use of the VP 
algorithm, with the initial EIR guess shown in the right panel of subfigure (b). The zoomed-in 
image corresponds to the ROI of the dashed rectangle. The grayscale window was [0,0.06]. The 
right plot shows the recovered and initial EIR. The SIR was accounted for in both cases. 
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(b) A = 1.0 X 10-^ A = 1.0 X 10-^, A = 1.0 X 10“^ a = 1000 


n 0.16 


10.08 





0.16 


0.08 


(c) A = 1.0 X 10“®, A = 1.0 X 10“^, A = 1.0 X 10-^, a = 10000 



0.16 


0.08 



Time (|i s) 



(d) A = 1.0 X 10-®, A = 1.0 X 10-^, A = 1.0 X lO"®, a = 10000 


Fig. 18: (a) and (c) Images reconstructed using the VP algorithm and different regularization 
parameters A with the same initial guess of EIR. a = 1000 for all reconstructions in (a), and 
a = 10000 for (c). The zoomed-in image corresponds to the ROI of the dashed rectangle. The 
grayscale windows were [0,0.16]. (b) and (d) show the corresponding recovered EIR. The SIR 
was accounted for in all cases. 
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(a) Images reconstructed by use of the conventional iterative method with different speeds of sound. The speeds 
of sound employed in the first, second, and third images are c = 1.542, c = 1.54, and c = 1.535, respectively. 
A = 1.0 X 10“^ for all reconstructions. The grayscale windows were [0,0.24]. 



(b) Images reconstructed by use of the VP algorithm with different speeds of sound. The speeds of sound employed 
in the first, second, and third images are c = 1.542, c = 1.54, and c = 1.535, respectively. A = 1.0 x 10“^, and 
a = 5000 for all reconstructions. The grayscale windows were [0,0.24]. 



(c) Recovered EIR 


Fig. 19: Reconstructed images and EIRs using different speeds of sound. The zoomed-in image 
corresponds to the ROI of the dashed rectangle. 
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